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ABSTRACT 

We combine spectra of VY CMa obtained with the short- and long- wavelength 
spectrometers, SWS and LWS, on the Infrared Space Observatory to provide a 
first detailed continuum spectrum of this highly luminous star. The circumstellar 
dust cloud through which the star is observed is partially self-absorbing, which 
makes for complex computational modeling. We review previous work and com- 
ment on the range of uncertainties about the physical traits and mineralogical 
composition of the modeled disk. We show that these uncertainties significantly 
affect the modeling of the outflow and the estimated mass loss. In particular, we 
demonstrate that a variety of quite diverse models can produce good fits to the 
observed spectrum. If the outflow is steady, and the radiative repulsion on the 
dust cloud dominates the star's gravitational attraction, we show that the total 
dust mass-loss rate is ~ 4 x lO~^M0yr~^, assuming that the star is at a distance 
of 1.5 kpc. Several indications, however, suggest that the outflow from the star 
may be spasmodic. We discuss this and other problems facing the construction of 
a physically coherent model of the dust cloud and a realistic mass-loss analysis. 

Subject headings: circumstellar matter - infrared: stars - stars: abundances - 
supergiants - molecular processes - stars: individual (VY Canis Majoris) 

1. Introduction 

VY CMa is a highly luminous, variable M2/3 II supergiant. The star's spectrum is 
dominated by a huge infrared excess. Le Sidaner & Le Bertrc (1996) have shown that the 
star is obscured, with Aj = 3.09 - 3.35, corresponding to reddening E(B-V) ~ 4.6. Only 1 
percent of the total luminosity is detected at optical wavelengths. The star's luminosity 
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has been estimated by a number of groups (cf Lada & Reid, 1978, Le Sidaner & Le Bertre, 
1996, Richards et al, 1998 ), but remains somewhat uncertain because its distance is not 
well-estabhshed. The consensus seems to converge on a distance of 1500 pc and a luminosity 
L = 5 X IO^Lq. 

The star exhibits strong outflows. Buhl et al. (1975) detected thermal emission in 
the vibrational ground state (v = 0, J = 2 - 1) 86.84686 GHz SiO line, which Reid and 
Dickinson (1976) interpreted in terms of a stellar radial velocity of 17.6 ± 1.5kms~^ with 
respect to the local standard of rest, and an expansion velocity of 36.7±2.0 kms~^. Neufeld 
et al. (1999) studied three thermally emitted mid-infrared water vapor emission lines in a 
high-resolution spectrum obtained with the Short Wavelength Spectrometer (SWS) on the 
Infrared Space Observatory (ISO), and found a mean radial velocity of order 20 ± 2kms~^ 
and a 25kms~^ outflow velocity. This outflow velocity is significantly lower than the 
velocities cited by Reid & Dickinson (1976) or the 32kms~^ that Reid & Muhleman (1978) 
reported for 1612 MHz OH maser outflow. 

Le Sidaner & Le Bertre (1996), who adopted a radial velocity of 20.5 km s~^, a 
galactocentric distance to the Sun of 8.5 kpc, and a rotational velocity about the Galactic 
center of 220kms~^ derived a distance to VY CMa, D = 2100 pc. Earlier work by Lada 
& Reid (1978) associates VY CMa with a molecular cloud complex, one edge of which 
appears to be illuminated by the star. Since the radial velocities of the complex and star 
are rather similar, the star and cloud appear to lie at the same distance, which Lada & Reid 
estimated at 1500 pc. On this basis they derived the star's luminosity as L = 5 x IO^Lq. 
The difference in the distance estimates of Lada & Reid on one hand, and Le Sidaner & Le 
Bertre on the other, appears to be based on somewhat different modeling of the Galaxy, 
since both estimates depend on roughly the same radial velocities. 

An independent means of estimating the star's distance is offered by the radial motions 
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of circumstellar masers. VY CMa exhibits intense OH, H2O and SiO maser emission 
indicating rather complex spatial distributions and outflow velocities, particularly in the 
H2O maser lines (Bowers, Claussen, & Johnston, 1993). Radio interferometric H2O maser 
observations that Richards et al. (1998) undertook after the pioneering observations by 
Bowers, Claussen and Johnston indicate substantial proper motions in the nine years from 
1985 to 1994, generally directed away from the assumed position of the central star. If the 
star lies at a distance of 1500 pc, these proper motions translate into velocities of order 8 
km s~^ at a distance of 75 mas from the star and 32 km s"^ at 360 mas. Such velocities are 
comparable to the radial velocities of the H2O masers, which exhibit a maximum expansion 
velocity of 35 kms~^ at 400 mas from the central position. 

The effective temperature for VY CMa appears to be Teg = 2, 800 K, agreed on, among 
others, by Monnier et al. (1999a), who adopt a distance of 1.5 kpc, and Le Sidaner & Le 
Betre (1996), who preferred a distance of 2.1 kpc to the star. The combination of effective 
temperature and a luminosity of 5 x lO^L© then fix the star's radius as i?* ~ 14 AU, though 
Monnier et a/.(1999a) cite a lower value ~ 2 x lO^L©, and a correspondingly lower stellar 
radius of 8.8 AU on the basis of their near-infrared observations ranging out to ~ S/im. 

High spatial-resolution, near-infrared observations of VY CMa have become available 
in the last few years. Danchi et al. (1994) used heterodyne interferometry at 11 /im, at 
intervals that covered both maximum and minimum luminosity phases of the star. For 
these respective phases, they model their data in terms of photospheric radii of 9.5 and 11 
mas and inner radii for the dust shell of 50 and 40 mas. The outer radius of 2.5 arcsec, 
remains the same for both phases. A photospheric radius of 10 mas would correspond to 
a stellar radius of 15 AU for a star at a distance of 1500 pc, in rough agreement with the 
conclusions of Lada & Reid. 

Wittkowski et al. (1998) obtained measurements of the circumstellar envelope's angular 
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diameter at 0.8, 1.28 and 2.17//m, which showed a nonspherical disk-hke structure with 
respective dimensions 60 x 83, 80 x 166 and 138 x 205 mas. The star itself was not resolved 
in their observations. Monnier et al. (1999a) succeeded in mapping the rather irregular 
dust cloud surrounding the star, at wavelengths of 1.65, 2.26 and 3.08 //m. They established 
that the hne of sight optical depths to the star at the respective wavelengths are 1.86 ± 
0.42, 0.85 ± 0.20 and 0.44 ± 0.11. 

Efstathiou & Rowan- Robinson (1990) derived a model for VY CMa, which explains 
the object's major features. They consider the central star to be embedded in a massive 
disk of dust and gas flowing out from the star at constant velocity. The disk flares out from 
the star at angles ±45° above and below the disk's central plane. The inner edge of the 
disk at radial distance ri has dust at a temperature Ti = 1000 K. The disk stretches to a 
radial distance r2 from the star, with ri/r2 = 0.006. At angles > 45° off the plane, radiation 
emitted by the star skirts over and under the disk and escapes into space. Efstathiou and 
Rowan-Robinson beheve that the hne of sight from Earth is inclined at about 43° to the 
disk's orbital plane, so that we are looking almost tangentially through the outer layers 
of the disk. Their model accounts for the low levels of visible radiation as well as the 
relatively high reddening observed by Monnier et a/. (1999a). The optical light we receive 
from the star is partly due to scattering, and is largely attenuated through absorption by 
circumstellar dust. The radius ri in the model of Efstathiou & Rowan- Robinson should 
roughly correspond to the figure of 40 - 50 mas given by Danchi et al., or ~ 70 AU, while 
their value of r2 could correspond to either the value cited by Danchi et al. ~ 2.5 arcsec, or 
to the extent found by Monnier et a/.(1999a) at 1.25 /xm, where the cloud extends beyond 
4 arcsec along some directions. These observations respectively yield ri/r2 ~ 0.02 to 0.01, 
only shghtly higher than the value suggested by Efstathiou & Rowan-Robinson. 

It is possible that the work of Efstathiou & Rowan- Robinson could now be extended 
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using methods developed by Chiang and Goldreich (1997, 1999) for calculating the spectral 
energy distribution of disks surrounding T Tauri stars. While these stars are quite different 
from VY CMa, and the disk surrounding VY CMa involves outflow rather than infall, 
several of the insights developed by Chiang and Goldreich apply to both kinds of dust 
shells. 

Mass loss rates cited for VY CMa by different authors range from 1 x lO~^M0yr~^ 
(Jura & Kleinmann, 1990) to 3.1 x IO-^MqJj;-^ (Danchi et al, 1994), with other workers 
citing somewhat different numbers. All agree, however, that the mass loss rate is of this 
general order. 

The structure of the present paper reflects the many significant uncertainties remaining 
about VY CMa, in particular the wavelength-dependent optical depth of the circumstellar 
cloud, the consequent acceleration and spatial density distribution of the outflow from the 
star, and the chemical and mineralogical composition of the dust grains. In sections 2 and 
3 we describe, respectively, the observations we carried out and the spectrum obtained. 
Section 4 discusses the current state of modehng and the respective weaknesses of a broad 
range of models that may provide excellent fits to the data but remain unphysical. Section 5 
introduces an unabashedly unphysical toy model of VY CMa, primarily for the purposes of 
examining the possible mineralogical constituents that might be expected in a satisfactory 
model of the VY CMa outflow. Section 6 exphcitly lists the physical weaknesses of this 
model, but indicates why the chemical nature of the various constituents might nevertheless 
be correctly identified. Sections 7 and 8 provide two different grain mass-loss rates from 
the star. Sections 9 calls attention to a minor anomalous spectral feature at ~ 7 /im, for 
which we found no satisfactory explanation. A summarizing discussion and our conclusions 
follow, respectively, in sections 10 and 11. 
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2. Observations 

The observations on VY CMa reported here were obtained using both the short- and 
long- wavelength spectrometers, SWS and LWS, on the Infrared Space Observatory, ISO 
(Kessler et al. 1996). Because several hours of observing time were required to cover the 
entire spectral range and because of scheduling conflicts the data were obtained on different 
ISO orbits. The wavelength ranges from 2.95 to 5.30 and 12.0 to 19.5 /xm were scanned on 
13 October, 1997. All other portions of the spectrum were scanned on the orbits of 19, 21 
and 24 November, 1997. The observations of November 21 duplicated and extended the 
2.95 to 5.30 //m observations of October 13 and were in good agreement. The wavelength 
range from 12 to 19.5 /xm was not revisited during the November observing runs, but our 
combined spectrum shows no discontinuities, suggesting that the star's spectrum did not 
appreciably change during this period. This supposition is given weight by the observations 
of Forrest, McCarthy and Houck (1979) who observed VY CMa from the NASA Lear Jet. 
They found the 20 /xm emission of VY CMa to be essentially constant in measurements 
carried out in January and November, 1976 and September, 1977. In particular, the last 
two measurements, taken ten months apart, differed by only 4%. Our ISO 20 //m flux of 
8 X 10~^^ Wcm~^/^m~^ is indistinguishable from their cited September 1977 value, taken 
twenty years earlier. This does not mean that the flux from VY CMa is steady. Monnier, 
Geballe & Danchi (1999b) report 8 - 13 //m variations on the order of ±20% over several 
years, but significant variations are not expected on a time scale of 6 weeks for this irregular 
or semiregular variable whose approximate period is 2200 days. 

The short wavelength spectrometer uses different apertures in scanning different 
parts of the spectrum. In the wavelength range from 2.38 to 12 /xm, the instrument views 
a 14" X 20" rectangular beam on the sky. At 12 to 29.5 jira the corresponding beam 
dimensions are 14" x 27", and at 29.5 to 45.2 //m, the beam size increases to 20" x 33". 
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These beams all were centered on coordinates a = 7''22"*58."3, b = -25°46'03" . The 
spectral resolving power A/5 A ranged from ~ 1000 to ~ 2000 in different portions of the 
spectrum. The total observing time on the star was 10,185 s. 

A group led by Michael Barlow observed VY CMa with the long wavelength 
spectrometer (LWS) on the night of 20 November, 1997. Their spectrum obtained with 
a circular ~ 80" diameter field of view covered the wavelength range from 45 to 200 //m 
in a run lasting 1930s (Barlow 1999). We made use of these data obtained from the ISO 
archives. 

Despite the difference in the sizes of the fields of view, we believe that even the smallest 
of the apertures used at the shortest wavelengths will have captured most of the radiation 
from the dust cloud enveloping the star. Our smallest field of view extends well beyond the 
faintest near-infrared regions mapped by Monnier el al. (1999a). We do not expect even 
the coldest portions of the dust cloud radiating in the mid- and far-infrared, to lie outside 
the fields of view of the ISO spectrometers. 

The flux integrated over our entire spectrum corresponds to ~ 6 x lO^^^Wcm"^. 
This may be compared to a total flux of ~ 7.1 x 10~^^ Wcm~^ expected from a star with 
bolometric luminosity 5 x IO^Lq at a distance of 1.5 kpc. We, therefore, adopt an infrared 
dust cloud luminosity of ~ 4 x IO^Lq. In the model of Efstathiou and Rowan- Robinson 
the outflow around the star occupies a disk that extends ~ 45° above and below its central 
plane. Such a disk occupies a field of view of ~ 9 sr as seen from the star. If the star's 
luminosity is 5 x IO^Lq, roughly three-quarters of the star's radiation will be absorbed 
by the disk and re-emitted isotropically in the infrared, while the remainder will escape 
unattenuated in directions surrounding the disk axis of symmetry. The total observed 
infrared flux is consistent with this model and somewhat lower than expected if the star 
were fully encased in a dust shell. 



-9- 



3. Results 

The reduction of the SWS grating spectrum, which ranges from 2.38 to 45 //m, was 
carried out with the SWS Interactive Analysis system using cahbration files from version 
7.0 of the ISO pipeline software. The LWS spectrum obtained from the ISO archives had 
been reduced with LWS pipeline 7. Figure 1 shows a composite of all available data. 

Figure 2. is a plot of the short-wavelength data out to 45 /xm, which shows a number of 
highly broadened features and also registers some of the strongest spectral hues. The hne 
spectrum between 29.5 and 45 iim has already been discussed by Neufeld et al., (1999). 

4. Weaknesses of Current Models 

Some of the first models of VY CMa were primarily descriptive. Early spectra of VY 
CMa already revealed strong silicate emission bands at 10 and 20 fim.. These motivated Le 
Sidaner & Le Bertre (1996) to use the "dirty silicate" model of Jones & Merrill (1976) in 
discussing their observations. 

The first computational model to deal with both the mineralogy and geometry of the 
circumstellar disk was due to Efstathiou & Rowan- Robinson (1990), who took optical depth 
effects into account and obtained a largely satisfactory fit to the spectrum of VY CMa. 
These models, though by now widely accepted in general outline, nevertheless entailed 
many simplifications and somewhat arbitrary parameters. One simplification was a ET^ 
dependence of the density of the disk on radial distance from the star. The second was 
the density dependence on the altitude angle, i.e. the angle above or below the central 
symmetry plane that the disk subtends at the star. Efstathiou and Rowan-Robinson 
assumed this to follow a power law up to some maximum "opening angle" of the disk, 0, 
though they note that this choice was "not based on any physical grounds but ... introduced 
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in order to investigate its effects on the observable properties of the system." The third 
was the assumption of spherically symmetric scattering by the grains. With these three 
caveats in mind, Efstathiou and Rowan- Robinson then computed a whole range of models 
with different opening angles as seen by an observer placed at different viewing angles - the 
choice of viewing angle being a fourth free parameter. 

By examining a large number of power-law fits and viewing angles, in this class of 
models, the best model fit, obtained on the further assumption that only "dirty silicate" 
grains were involved, suggested a central star surrounded by a disk with opening angle 
45° above and below the plane, with the terrestrial observer viewing the star at an almost 
glancing angle of 43° above the disk, so that the line of sight to the observer still cuts 
through an appreciable amount of obscuration. Radiation the star emits at angles greater 
than 45° above or below the disk propagates unattenuated into circumstellar space. 

An important limitation of the Efstathiou and Rowan- Robinson model is that it 
neglects the coupling between opacity and dynamics. The dust grains absorb the star's 
radiation most strongly at small distances from the star, where the outflow accelerates and 
the density correspondingly drops more rapidly than R~^. Only at large distances from 
the star, where the radiation pressure is low does the outflow reach a constant terminal 
velocity. But at these distances the irradiation also has dropped, and the grains evidently 
will hardly radiate at all. The assumption of constant velocity, therefore, may reflect a 
quite unphysical simplification. 

That Efstathiou and Rowan- Robinson, nevertheless obtain a good fit to the observed 
spectrum may be seen as a sign of ambiguity. A good model fit to the spectrum does not 
necessarily imply a physically correct model. 

A first attempt to reduce this ambiguity by providing a model which couples radiative 
transfer to dynamics was published by Ivezic & Elitzur (1995). These authors take 
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into account not only the radiation pressure acting on the dust, but also the countering 
gravitational attraction of the star. Zubko and Elitzur (2000) have recently applied this 
model to the outflow of water and dust from the star W Hydrae. Their model, however, is 
spherically symmetric and would have to be modified to deal with the disk-shaped outfiow 
now widely assumed to characterize VY CMa, or with the plume of material that Monnier 
et al. (1999a) suggest may be flowing out of the star, spraying its surroundings like a 
rotating garden hose. 

The present authors are not aware of any computer model currently available, that has 
sufficient generality to permit a full description of the star. 

5. A Toy Model Outflow 

In the present paper, we resign ourselves to not being able to provide a model that 
convincingly couples radiative transfer to momentum transfer in a realistic outfiow of 
condensing grains, and limit ourselves to the presentation of our improved spectrum and 
the pursuit of a somewhat complementary descriptive approach. We set up a toy model 
that has many of the same physical shortcomings as others but may help to illuminate the 
mineralogy somewhat better, given the much more detailed spectral energy distribution 
now available. In doing this, we set aside the inspired geometric approach of Efstathiou and 
Rowan- Robinson (1990), and the elegant dynamic approach of Ivezic and Ehtzur (1995), in 
what is clearly no more than an interim sketch. 

Our toy model is spherically symmetric and optically thin to radiation re-emitted 
by the grains. Both these assumptions are easily challenged, but even these unrealistic 
assumptions can yield a remarkably good fit to the spectrum, and provide a number of 
insights to be enumerated below. The assumptions made permit us to use a straightforward 
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version of the radiative transfer code MODUST, which incorporates information on the 
mineralogical composition of the dust and permits the introduction of several different 
mineralogical components. 

In Figures 1 and 2 we present the MODUST fit obtained for the spectrum of VY 
CMa. To compare results to other workers in the field, we have based all of these models 
on an assumed stellar luminosity of 5 x lO^L© and a distance to VY CMa of 1.5 kpc, 
even though the flux we observe corresponds to a luminosity roughly 20% lower, for a 
star at this distance. This is consistent with our assumption that, as in the disk model of 
Efstathiou and Rowan- Robinson, 20% of the flux escapes into space unattenuated, while 
80% is absorbed and re-emitted in the infrared for an infrared luminosity of 4 x IO^Lq. 

Our model assumes a i?~^-dependence of the number-density of grains on radial 
distance from the star. It invokes two major mineralogical constituents, but also illustrates 
the lesser contributions due to a potential admixture of two minor components. The major 
constituents are metallic iron and amorphous silicates, the minor components are crystalline 
sihcates. We included iron mainly because sihcates are poor absorbers at short wavelengths 
where VY CMa exhibits strong absorption. The minor components show that crystaUine 
silicates can at best provide only small contributions to the overall spectrum. 

Our choice of constituents is motivated by condensation sequences applicable to an 
oxygen-rich star. If we assume that the dust is formed from a vapor of solar composition 
during the early stages of the outflow, metallic iron flrst condenses at T ~ 1460 K, 
followed by calcium-rich pyroxenes and nearly pure magnesium olivines at T ~ 1450 K 
(see also Whittet 1992). Nuth et al. (2000) argue that, in a thermodynamic equilibrium 
condensation sequence, magnesium silicates will condense well before iron silicates become 
stable at around 600 K. This makes it highly unlikely that Fe-containing silicates will 
form on time scales as short as those considered here, since the solid-state reaction rates 
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between Fe-metal and magnesium silicate grains are slow. Even when a Mg-Fe-SiO-vapor is 
condensed, the presence of Fe-atoms in the resulting grains excludes the presence of Mg and 
vice- versa (Rietmeijer et ai, 1999). Given these findings, we assume that the amorphous 
silicates are pure Mg-rich olivines. 

Our best fit to the data is shown in Figs. 1 and 2. The grain parameters used in the 
model are given in Table 1. We have used the optical properties of amorphous silicates 
that were available from laboratory research, and emphasize that most of the parameters 
were chosen primarily to obtain the best possible fit. The temperature of the metallic iron 
particles in our model ranges from 630 to 110 K, dropping from the inner to the outer 
shell, corresponding to distances from the star indicated in the second column of Table 
1. In general, one might expect the formation of FeO and FeS instead of Fe at the cited 
temperatures. However, oxidation of the metal grains, which would give rise to a strong 
emission peak at 23 /im, appears to be forestalled by an earlier incorporation of most of the 
free oxygen in H2O and other molecules. Similarly, the presence of triohte (FeS) would give 
rise to several strong peaks in the near- and mid-infrared, which are not discerned in the 
spectrum. 

Compared to iron, silicate grains absorb starlight less efficiently and re-radiate more 
efficiently at long wavelengths. They can therefore condense closer in to the star and still 
have a lower temperature. As Figure 2 demonstrates, the potential contributions from 
crystalline forsterite (Mg2Si04) and crystalline enstatite (MgSiOs) are quite minor, and it 
is not clear that these components are actually present. We have included their spectra 
in our model primarily to indicate the wavelengths at which contributions from these and 
similar dust species would be significant and the extent to which small admixtures of 
different grain materials could affect the spectral fit. Because the mid-infrared spectrum of 
VY CMa exhibits small deviations on top of a smoother continuum emission, the presence 
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of crystalline silicate dust is not unlikely; however, we cannot identify the actual carriers. 
In this respect the dust disk around VY CMa is quite different from the disk surrounding 
the young star HD 100546 which exhibits strong forsterite emission (Waelkens et al, 1996). 

Two limitations now need to be mentioned: 

1. If we adopt the Solar System elemental abundances of Anders & Grevesse (1989), 
the number densities of iron, magnesium and silicon atoms should all be very nearly 
identical. The atomic weight of iron is 56 atomic mass units, amu. Rietmeijer et al, (1999) 
find that condensates from Fe-Mg-SiO-H2-02 vapor form MgO-Si02 sohds that can take 
on the form of pure amorphous silica Si02 and also pure MgO, but tend to cluster around 
compounds like MgSiOs, with equal numbers of magnesium and silicon atoms. For MgSiOa, 
the molecular weight is 100 amu, meaning that there is one atom of silicon plus one of 
magnesium per 100 amu of silicates. The iron-to-sihcate mass ratios of 55 to 45 % that we 
find in the best-fitting Toy Model with parameters listed in Table 1, would then lead to 
elemental ratios of Fe:Si:Mg = 2:1:1. At face value, the iron abundance, therefore, appears 
to be rather higher than expected though not too surprising, given the crudeness of the toy 
model. 

The importance of including iron in our models needs to be stressed. Iron grains 
should naturally condense at high temperatures and absorb well at near infrared frequencies 
where the star emits much of its radiation but pure silicates absorb poorly (cf. Jones and 
Merrill, 1976). The choice by Efstathiou and Rowan- Robinson to base their model on "dirty 
silicates" may lead to similar spectra, but since the mineralogical nature of the "dirt" in 
these silicates comes from empirical fits found by Jones and Merrill, it does not readily 
lend itself to calculations of mass loss rates, or an understanding of the chemistry of the 
outflowing material. 

2. A somewhat suspicious feature of our model is the sharp rise in the spectrum at 
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short wavelengths. Out to about ~ 5 //m it mimics the Wien distribution of a blackbody 
at temperature ~ 400 K. This suggests that VY CMa is hidden from direct view by a 
dust photosphere, which need not surround the star entirely but at least prevails for 
rays passing radially outward through the disk. The existence of such a photosphere is 
consonant with the high near- infrared optical depths measured by Monnier et al. (1999a), 
which indicate that most of the starlight is absorbed by dust along the line of sight. If 
the dust-photospheric portion of the spectrum is close to blackbody radiation, then the 
chemical properties of the radiating dust become less important, and the sharp rise in the 
observed spectrum at short wavelengths no longer needs to be attributed to any particular 
constituent. However, iron grains do provide a means for readily obtaining substantial 
optical depths at wavelengths emitted by the star. 

We examine this question in the next section. 

6. Weeiknesses of the Toy Model 

The similarity of the iron emission curve in Fig. 2 to blackbody emission indicates that 
the optical depth of the iron content of the model needs to be re-examined. Our model 
assumed spherical iron grains of radius a, with a size distribution n(a) = n(ao)(ao/a)^'^ and 
a radial distribution proportional to R~^i where R is the distance from the star. This radial 
distribution corresponds to radial outflow of a homogeneous cloud at constant velocity, and 
implies that shells of equal thickness at any distance from the star contain equal mass. 
Calling the inner radius of the iron dust cloud Rq, we obtain a particle number- density 
distribution function 

n(a, R) = n(ao, RQ){a^/af-\R^/Rf (1) 



-16- 



This leads to a total iron dust mass 
M 

or 
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M = —n{ao, I^)al 'pRl{alll - aJ^J{Rmax - Ro)^ (3) 

Here p is tlie density of individual iron grains, and Q is the solid angle the dust cloud 
subtends at the star. 

The fraction of the radiation F absorbed by these grains, i.e. the fraction of the star's 
light intercepted in a radial transit through the cloud for the optically thin case is 
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In the present case, F can exceed unity at some wavelengths, indicating that optical depth 
effects need to be taken into account. We can rewrite this equation as 
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where the Mie theory value for the absorption efficiency can be taken as Qabs — 
87ra/A/m[(m^ — l)/{rn? + 2)], A is the wavelength, and Im indicates taking the imaginary 
part of the function in square brackets for a complex index of refraction m. The ratio of 
absorption to mass is then 



M R^RjnaxP^^ 



Im 



m^ — 1 
m2 + 2J 



(6) 



For the Toy model, we used p ~ Sgcm""*, Rq = 1.8 x lO^i?© = 1.26 x lO^^'cm, 
Rmax = 4.5 X IO^Rq = 3.15 X 10^^ cm. The total iron dust-mass in the model is 
M = 9.4 X ~ 1.9 X 10^^ g. Using the Efstathiou & Rowan- Robinson (1990) 

configuration, the dust disk occupies a field of view Q ~ 9 sr as seen from the star. This 
leads to 



12.5 
IT 



Im 



-m - 1" 
.m^ + 2. 



(7) 



-17- 



where the wavelength A is measured in micrometers. Only ~ 7% of the star's radiation 
comes from wavelengths as short as 0.73 //m, where F has as high a value as ~ 2.5, 
indicating virtually total absorption at these wavelengths in a thin shell closest to the 
star. At l//m, the complex index of refraction is m ~ 3 + 4i (Pollack et al., 1994), so that 
the imaginary term is ~ 0.1. At wavelengths longer than ~ 1 //m, i.e. for roughly three 
quarters of the energy output of a star at temperature 2800 K, the imaginary term rapidly 
diminishes to values less than 0.1, and F even more rapidly diminishes to a value below 
1.25. By A ~ 1.5 //m, roughly the point beyond which the star still emits half its light, F 
has dropped to ~ 0.5 and F should become a fairly good indicator of the fraction of the 
radiation absorbed. Taken together, these figures indicate that virtually all the starlight 
at wavelengths short of 1 //m, and more than half the light between 1 and 1.5 //m is fully 
absorbed in this model. 

Since the number of grains in shells of equal thickness around the star remains constant, 
and the outflow velocity is taken to be constant, the differential absorption F{R) dechnes 
with increasing distance from the star as R^^. This has two immediate consequences for 
a cloud that stretches radially from a distance Rq out to much larger distances — in our 
model out to 25 Rq. Three quarters of whatever energy the cloud absorbs, is already 
absorbed by the time radiation reaches a distance 2i?o, so that the inner parts of the cloud 
closest to the star act as though there were a relatively thin dust photosphere. And since 
the absorption integrated over most of the star's spectral energy distribution function over 
this radial distance interval from the star is less than unity, the MODUST program will at 
least give roughly reasonable results. 

The upshot of this argument is that the type of dust-photosphere which emerges from 
the model easily explains the steep blackbody-like rise in the VY CMa spectrum observed 
around 4 /xm, while also providing a superior fit at long wavelengths where a blackbody 
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photosphere would clearly give too low a flux. 

For the spectrum due to amorphous silicates, Draine & Lee (1984) plot Qabs values, 
which indicate a range of /m[(m^ — l)/(m^ + 2) < 0.04 for all grains smaller than 1 at 
essentially all wavelengths spanning the star's spectral energy distribution. For the assumed 
silicate grain density of 3 g cm~^ in our model, this implies that 



^ 22.5 , 
F = —-7m 
A 



<1 (8) 



.m2 + 2. 

so that the absorption in the region between Rq and 2i?0i where most of the radiation is 
absorbed is always lower than 0.5 for grains 1 /im in diameter and smaller. The MODUST 
program, therefore, again can be expected to yield reasonable values. 

One might, however, still wonder why the Toy model silicate emission should so closely 
resemble that of the observations. The answer to this is that silicates absorb roughly 
equally well at 1.5 /xm, where the star's spectral energy distribution is centered, as at 
10 /xm, where the blackbody photosphere or the iron photosphere of our Toy model have 
their respective spectral peaks. The curves drawn by Draine & Lee (1984) show that 
Im[{m^ — 1)/ {rri^ + 2)] rises to a value of ~ 0.5 around 10 //m and is remarkably independent 
of grain size for the range a = 0.01 to 1 fim. This means that the fraction of the radiation 
absorbed F oc ImKni^ — l)/(m^ + 2)]/A, and the implied grain absorption, is roughly the 
same for photospheric emission from a star at temperature 2800 K as for emission from a 
dust photosphere of comparable luminosity. The MODUST program, therefore, fortuitously 
provides roughly correct results either way. 



7. Toy Model Mass-Loss Estimate 



Despite the limitations just enumerated, our model should be able to yield a rough 
estimate of the dust-mass loss from VY CMa. The total dust mass in the Toy model is of the 
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order of 1.6 x lO^^M©. The model reaches out to a distance of ~ 1550i?* = 3.3 x 10^^ cm, 
rather larger than the optical and near-infrared extent of the cloud measured by Monnier 
et al. (1999a), who placed the outer edge of their observed dust disk at around 4 arsec, or 
~ 10^^ cm. For an outflow velocity of 25-30 km s~^ our model disk would have to replenish 
itself every 3500 — 4200 yr. The minimum dust- mass-loss rate thus would become about 
~ 3.8 — 4.6 X 10~^Mq yr^^ if we assume that all of the disk replenishes itself this rapidly. If 
the gas-to-dust ratio in the outflow is of order 100, the grain mass loss we derive translates 
into an overall mass of loss of M ~ 4 x lO^^Mgyr"^. Both values are a factor of roughly 4 
higher than estimates of the mass-loss obtained by other observers. 

Sopka et al. (1985) derived a grain mass loss of ~ 8 x lO~^M0yr~^. This estimate 
was based on two considerations: an assumed dust opacity law that ranges all the way 
from the ultraviolet to submillimeter wavelengths and a 400 /xm flux of 10 Jy - which can 
be compared to our best estimate of 15 Jy. The remaining difference in our estimates 
therefore amounts to a factor of order 3, which is understandable in view of the appreciable 
differences in approaches and assumptions. 

Danchi et al. (1994) estimated a gas mass loss of 3.1 x lO^^Mgyr"^, though their 
results cannot be directly compared, because they assumed a mixture of silicates and 
graphite, and a temperature at the inner radius of the dust shell as high as 1536 K. The 
dust temperature they use in interpreting their visibility curves appears to give too high 
a temperature and would result in higher near-infrared emission than we detect. They, 
further, assumed a gas-to-dust mass-ratio of ~ 200, meaning that their dust mass loss 
would amount to only ~ 1.5 x 10~^MQyr~^. 

Our Toy model estimated mass loss is also rather higher than the mass loss of 
1.2 X lO~^M0yr~^ that Netzer & Knapp (1987) inferred for an adopted distance of 1900 pc, 
using CO (J = 2 - 1) data gathered by Zuckerman & Dyck (1986). Since this value is based 
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on CO, rather than dust outflow, a direct comparison is also rather difficult. 

8. An Alternative Approach to a Mass-Loss Estimate 

A totally independent mass-loss estimate, rather simpler than those based on the more 
detailed models described above, can be obtained on two assumptions. The first is that all 
of the observed infrared emission is emitted by dust and gas flowing out of the star with 
identical outflow velocities. The second is that the outward directed radiation pressure on 
the gas and dust appreciably exceeds the gravitational attraction to the star. 

These assumptions would not apply to re-emission from dust in a gravitationally 
bound photosphere. The assumptions then specifically refer to a steady-state outflow as in 
our Toy model, but might not hold if part of an emitting dust photosphere could remain 
gravitationally bound and possibly re-contract in a variable star like VY CMa. 

To calculate the mass loss under the stated assumptions, we let the infrared luminosity 
of the dust cloud be Ljr. Then the rate at which outward directed momentum is deposited 
by the star into the circumstellar disk is Ljr/c. Since the terminal velocity of the gas is 
known to be of order v^o = 25kms^^, the total mass accelerated from some initial velocity 
< Voo per unit time amounts to 

M = , , ~ ^ ^ , ~ 3 X 10-^Meyr-^ (9) 

c(voo - vo) 7.5 X 1016 (cm/s)' ^ ^ 

Under the stated assumptions the numerical value on the right is a minimum mass 
loss. If there were a substantial injection velocity vo, the mass loss would be higher. 

A more complete analysis that takes gravitational attraction into account leads to 
a more complex conclusion. Let the acceleration of the cloud take place between radial 
distances R and R + AR from the star, and the outfiow injection velocity at R be vq. The 



- 21 - 



acceleration then is 



MR^^-^^ (10) 



where is the mass of the star, and M the mass of the accelerating cloud. Let v = vq^ t^v 
be the velocity reached at distance R + Ai?. We readily integrate this with respect to time. 
For R ^ AR, we obtain 

M _ LiR M^MG 

Ai?/[^;o + (A^;/2)] cAv R^Av ^ ' 

The denominator on the left side of the equation is roughly the time during which the 

cloud is accelerated to velocity vq + Av, so that the full expression on the left gives the 

approximate mass loss rate. The first term on the right corresponds to the expression on 

the right of equation (9). If the gravitational attraction given by the second term on the 

right is substantial, the mass loss rate not surprisingly decreases. (For a rigorous discussion 

of this topic, see Ivezic &Elizur, 1995.) 

Equation (11) tells us that the very simple estimate of outflow given by equation (9) 
can only yield an upper limit to the mass loss. For VY CMa this upper limit, however, 
appears to be only a factor of 2 to 4 higher than values based on the arguments of Sopka et 
al. (1985), Danchi et al. (1994), or Netzer & Knapp, (1987). 



9. An Anomalous Spectral Feature at 6.7 to 7.5 /xm 

While we have been able to model most of the continuum spectrum of VY CMa 
through emission by warm dust, an anomalous feature at 6.7 to 7.5 fim has defied ready 
explanation. We checked whether water vapor fluorescence might contribute, but find that 
it most probably does not. 

Crovisier & Encrenaz (1983) and Crovisier (1984) investigated the possibility of H2O 
fluorescence in comets, and focused attention on a few optically-excited transitions that 
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appear of special significance. Gonzalez- Alfonso & Cernicharo (1999a) considered such 
processes also for circumstellar envelopes. 

Radiative absorption can raise a water molecule from the ground state, = 000, 

to the states 001 and 010 which, for illumination by cool stars, dominate by approximately 
an order of magnitude. Other transitions are less favored because the star emits too few 
photons at the higher required excitation frequencies or because the molecular Einstein B 
coefficients are too low. 

The sequence of events that could lead to the emission of significant radiation at ~ 7 
microns involves an initial excitation 000 001, and subsequent decay 001 010 followed 
by 010 000. If this process were in effect, we would expect to also detect an absorption 
feature due to the 000— s- 001 transition at 2.66 /im and an emission feature at ~ 4.36 /im 
due to the 001— >^ 010 transition. Neither of these is apparent. 

A second effect described by Gonzalez- Alfonso & Cernicharo (1999a,b) involves solely 
excitation from the ground vibrational state 000 to 010 and return to the ground state. 
In the downward transition from 010, some of the molecules end up in relatively excited 
rotational states, from which they subsequently decay through purely rotational emission. 
In the downward rovibrational transition to a higher rotational state the energy of the 
emitted photons is reduced and their radiation is shifted toward longer wavelengths. A 
typical spectrum will therefore show an absorption dip shortward of 6 /im, and an increase 
beyond 6//m. Gonzalez- Alfonso and Cernicharo (1999b) have modeled this effect for 
evolved stars, but their computations show that it is always more modest and restricted to 
individual spectral lines rather than to a superposition of faint lines so dense as to produce 
a strong continuum. In addition, the energy in the hump observed at 7 /^m far exceeds any 
absorption below 6 /im, and also stretches to appreciably longer wavelengths than the model 
of Gonzalez- Alfonso and Cernicharo predicts. We therefore conclude that this feature is not 
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produced by water vapor. 

Tsuji et al. (1998) have tentatively identified a broad emission feature tliey observed in 
tlie M2III giant (3 Peg at ~ 7.6 //m witfi SiO. Tfiis molecule should also exhibit a feature at 
4.0 lira.. We detect no corresponding structure at 4 /xm, and our 7 //m feature falls at shorter 
wavelengths. 

Dust species such as the silicate diopside (e.g. CaMgSi206, Koike & Shibai 1998) and 
crystalline spinel (MgAl204, Chiwara et al. 2000) exhibit a broad band at 7 yum (for spinel, 
this is not the case for a synthetic sample, but is true in natural white and pink spinel 
found in Burma and Sri Lanka). This band is accompanied by stronger bands at longer 
wavelengths. Amorphous diopside, which also exhibits the broad silicate bands at 10 and 18 
//m due to Si-0 stretching and 0-Si-O bending, is a particularly promising candidate. More 
detailed future laboratory analysis will permit testing these minerals as possible carriers of 
the 7 //m band. 

To date, we find ourselves unable to explain the 7 /xm feature, though we have no doubt 
that it is real. The mismatches at longer wavelengths, e.g. at 25 and 34 /xm can result 
from a slightly different chemical composition, ellipsoidal shape, or temperature of potential 
crystalline silicate components. 

10. Discussion 

In fitting the spectrum of VY CMa and estimating the star's mass loss, we have 
assumed that the radiating circumstellar dust cloud is fairly homogeneous. However, the 
patchy emission mapped by Monnier et al. (1999a) suggests the existence of a plume of 
material possibly ejected in different directions as the star rotates. 

Ejection of individual plumes or blobs is also suggested by the 22-GHz water maser 
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data of Richards et al. (1998) who found a hnear relationship between radial distance from 
the star and radial displacement over the nine year period from 1985 to 1994. This relation 
persists from 75 mas out to 350 mas (~ 100 to 500 AU) and covers the velocity range from 
~ 8 to ~ 32 kms~^. Remarkably, there is no apparent correlation between the observed 
line of sight velocity of the masers and their displacement from the star, but the linearity 
between proper motion and displacement extrapolates back to an intercept at zero distance 
from the star. Most workers in the field attribute the acceleration of both the thermal 
gas cloud and the masers to drag by dust grains propelled outward by radiation pressure 
from the star. Richards et al. find that the shell of water masers is elongated in the same 
direction as the infrared emission from the dust, and that the inner radius of maser emission 
appears to lie just outside the inner dust shell radius. The dust cloud and water masers, 
therefore, seem to be intimately related. However, the radiation pressure drops off as the 
inverse square of the distance, whereas a linear velocity /distance relation, r = ar, where a 
is a constant, requires an acceleration, r — ar — a^r, proportional to distance from the star. 
Even an optically thick cloud whose surface density declines as it expands radially outward 
from a star will not produce this high an acceleration as a function of distance. 

The linearity of maser velocity with displacement from the star can be most easily 
explained if all of the clouds producing maser emission were ejected at one and the same 
epoch and some were accelerated to higher velocities than others. Moreover, since the 
terminal velocity of the masers is comparable to the velocity spread Neufeld et al. (1999) 
observed in the thermal gas component, this would imply that the general stellar wind 
may similarly consists of individual blobs accelerated to a range of different velocities. 
This again implies that the bulk of the acceleration occurs where the radiation pressure 
is highest, near the inner radius of the dust cloud, in a dust shell that is optically thick 
at visual and near-infrared wavelengths but may vary in mass, i.e. surface density. The 
radiation pressure then accelerates low-mass blobs more rapidly and to higher terminal 
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velocities than higher-mass blobs. 

The imphcation of this view is that individual blobs ejected from the star are 
independently accelerated. This adds to the uncertainties still accompanying all models. 

If the lowest density blobs have optical thickness unity, r ~ 1, for photospheric radiation, 
and grains of density p ~ 4gcm~^ and radius a ~ 0.5 //m occur with a gas-to-dust ratio 
T^g/iTid ~ 100, the blobs' column densities a would be cr ~ rdprng/rrid ~ 4 x 10~^gcm~^, 
where d is the grain diameter. The radiation pressure at a distance R ~ 600 AU ~ 40i?* 
from the star is L^/At^B?c ~ 6 x 10~^dyn cm~^. This leads to an acceleration of 
~ 1.6 X 10~^cms~^, so that a velocity of 30kms~^ is reached in i ~ 2 x 10^ s, or ~ 60 
years. The distance the blob covers in that time is ~ 200 AU, which is sufficiently short 
that expansion and change in optical depth are quite minimal. The acceleration of the blob 
then stops primarily because material ejected from the star somewhat later begins to cool, 
form grains, absorb light, and thus cuts off the radiation from material at greater radial 
distance from the star. 

If a blob at distance R = 600 AU from the star, with column density cr ~ 4 x 10~^ gcm~^ 
is ejected at intervals of t ~ 2 x 10^ s and surrounds the star in solid angle of roughly 9 sr, as 
in the Efstathiou & Rowan- Robinson model, the mass outflow will amount to 5 x lO^^gyr"^ 
or ~ 3 X lO^^Mgyr"^, in agreement with our steady-state-outflow model estimate. 

Once a blob is accelerated to its terminal velocity the cloud material continues to drift 
outward, shaded by material that was ejected later. This shading is consistent with the 
observation by Monnier et al. (1999a) that the optical depth of the cloud is of order 2 at 
1.65 //m, and presumably higher at shorter wavelengths where most of the star's radiation 
is emitted. An implicit consequence of the argument presented in sections 6 and 8 is that 
infrared radiation comes only from the accelerating portions of the outflow. Portions of 
the cloud flowing outward at constant velocity do not experience radiative acceleration and 
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neither absorb starlight nor emit infrared radiation. 

11. Conclusions 

We have used a very simple model to obtain a rather good spectral fit to the infrared 
spectrum of VY CMa. This and the inferred mass loss may constitute a good first 
approximation despite the many uncertainties about this unusual star, and the many 
weaknesses of our model. 

Major constituents of the dust disk surrounding VY CMa appear to be amorphous 
silicates, probably with a significant contribution by iron. A good spectral match obtained 
to the near-infrared excess with metallic iron, Fe, may be fortuitous, though it would not 
be surprising, since iron particles should be first to form when a gas with solar abundance 
cools (see e.g. Whittet 1992). Prom condensation experiments by Nuth et al. (2000), it 
is apparent that amorphous magnesium silicates are also likely to condense in gases rich 
in silicon, oxygen, and magnesium. Small admixtures of crystalline enstatite and forsterite 
cannot be ruled out. At least a fraction of the grains appears to be several microns in 
diameter to account for the strong emission in the submillimeter and millimeter band. 

Our estimated mass-loss rate of order 4 x 1O~^M0 yr~^ is somewhat higher than values 
cited by other authors using rather different assumptions. However, all of the estimates 
agree that this stage of the star's evolution cannot last much longer than another ten to 
twenty thousand years, an interval during which 5 to 10 solar masses would be shed. Unless 
VY CMa is exceptionally massive, this is as much mass as the star may be expected to lose. 

As discussed in detail in section 4, even the most thoughtful models presented to 
date remain quite unphysical. To date, no model has been devised that meets most of 
the criteria, namely: (a) that it couples the dynamics to the radiative transfer, (b) that 
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it incorporates a realistic condensation sequence based on laboratory studies, (c) that it 
makes use of chemical abundances that are credible, and (d) that its geometry be plausible. 
Most models fail on all except perhaps one or two of these four criteria, and so we have to 
grope our way forward slowly. 

We have outlined the enormous amount of work that will be required before the nature 
of the circumstellar dust cloud around VY CMa is fully understood. Better computational 
tools, able to deal with a variety of chemically differentiated, optically thin and thick 
grain components and a more complex cloud geometry will be required and would lead 
to more realistic physical models. As input to such models, the condensation sequence of 
different chemical and mineralogical species in an oxygen-rich environment will need to 
be determined from laboratory studies. These could serve as input to computer models 
that couple an irradiated cloud's dynamical behavior to its density and temperature of 
condensed material in spherical, disk shaped or jet geometries. 
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Figure Captions 



Fig. 1- Overview of spectral data on VY CMa, obtained with ISO with the short wavelength 
spectrometer SWS and the long wavelength spectrometer LWS (solid jagged curve). Shown 
as triangles, are IRAS data at 12, 25, 60, and 100 //m, and ground-based submillimeter and 
millimeter data (Sopka et al. 1985; Marshall et al. 1992; van der Veen et al. 1995). Our 
estimate of the star's reddened flux is shown in the lower smooth curve, while the toy model 
best fit for stellar plus circumstellar emission is shown in the upper smooth curve (see text). 

Fig. 2 - Infrared Continuum Spectrum of VY CMa (0-50 jim) . The dash-dotted line gives our 
toy model fit (see text). The shaded regions provide the relative contributions by different 
dust components. The designation 'Cr.' stands for 'crystalline'. 
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Toy Model MODUST: Fit Parameters 



Parameters for the fit shown in Figure 1 and 2, together with their references 



Stellar parameters: 


• distance d = 1500 pc 








• radius R* = 3000 R© 








• effective temperature T^s = 2800 K 






• luminosity L* = 5 • 10^ 


Lo 


General dust parameters: 


• density law: N(R) = N^,,, (R/R*)-^ 






• radius law: n(a) = na,o 


a/ajjjjjj) 




MetaUic Fe 


Am. Sil. 


Cr. Mg2Si04 


Cr. MgSiOa 


Inner radius (R*) 


60 


25 


60 


20 


Outer radius (R*) 


1500 


1550 


1000 


500 


Nr,o (gCm-3) 


3 • 10-20 


9 • 10-20 


2.7 • 10-22 


4.5 • 10-20 


Particle radius, a, (//m) 


0.01 - 1 


0.01 - 5 


0.01 - 1 


0.01 - 1 


Grain density (gcm"^) 


7.87 


3.10 


3.33 


2.80 


Dust temperature (K) 


110 - 630 


50 - 440 


90 - 360 


60 - 170 


Total dust mass (M©) 


9.4 ■ 10-3 


6.6 ■ 10-3 


5.4 • 10-^ 


5 ■ 10-^ 


Relative percentage 


55 % 


41 % 


1 % 


3% 



Optical properties used: 



Amorphous silicates: Jagcr et al. (1994) & Mutschke et al. (1998) 
Crystalline sihcates: Jager et al. (1998) 
Metallic Fe: Pollack et al. (1994) 
Notes: 

'Cr.' indicates crystalline dust, while 'Am.' indicates amorphous dust. 
The internal grain density is not a free parameter but the density of the mineral 
studied in the laboratory to derive the optical constants of the respective species. 
All other parameters were chosen to provide a best fit to the data. 
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